Bullova et al. Figure 5 and Supplementary Figure 5
Jiacheng Zhu
2024-12-12
Load Data
Single-cell data
set.seed(1234)
# Load Data
WT_Mouse <- readRDS("./01.2_main_E17andPN_WT_AddMeta_Annotated.RDS")
WT_Mouse$NewCelltype <- ""
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 0)] <- "(0) E ADR Chromaffin"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 1)] <- "(1) PN NOR Chromaffin"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 2)] <- "(2) PN ADR Chromaffin"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 3)] <- "(3) Endothelial"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 4)] <- "(4) Macrophage"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 5)] <- "(5) E NOR Chromaffin"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 6)] <- "(6) E Glia|SCP"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 7)] <- "(7) Cortex ZG"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 8)] <- "(8) Cycling TH+"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 9)] <- "(9) Cortex X Zone"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 10)] <- "(10) E Neuroblast"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 11)] <- "(11) PN Sustentacular"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 12)] <- "(12) Cortex ZF"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 13)] <- "(13) Macrophage"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 14)] <- "(14) Cortex ZG"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 15)] <- "(15) Dendritic cell"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 16)] <- "(16) T cell"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 17)] <- "(17) NK cell"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 18)] <- "(18) Mesenchymal"
WT_Mouse$NewCelltype[which(WT_Mouse$seurat_clusters == 19)] <- "(19) Dendritic cell"
WT_Mouse$DetailedAge <- ""
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("E17nTAM_11_5", "E17nTAM_8_5"))] <- "Embryonic (E17)"
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("P5"))] <- "Postnatal (P5)"
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("P14"))] <- "Postnatal (P14)"
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("WT_AG_3m"))] <- "Postnatal (P90)"
WT_Mouse$DetailedAge[which(WT_Mouse$wenyu_names %in%c ("AdultAG"))] <- "Postnatal (Aged)"
WT_Mouse$DetailedAge <- factor(WT_Mouse$DetailedAge, levels = c("Embryonic (E17)", "Postnatal (P5)", "Postnatal (P14)", "Postnatal (P90)", "Postnatal (Aged)"))
WT_Mouse$GrossAge <- ""
WT_Mouse$GrossAge[which(WT_Mouse$wenyu_names %in%c ("E17nTAM_11_5", "E17nTAM_8_5"))] <- "Embryonic (E17)"
WT_Mouse$GrossAge[which(WT_Mouse$wenyu_names %in%c ("P5", "P14", "WT_AG_3m", "AdultAG"))] <- "Postnatal"
WT_Mouse$GrossAge <- factor(WT_Mouse$GrossAge, levels = c("Embryonic (E17)", "Postnatal"))
NewCelltypes <- sort(unique(WT_Mouse$NewCelltype))
NewCelltypesOrdered <- NewCelltypes[c(1,2,13:20,3:12)]
WT_Mouse$NewCelltype <- factor(WT_Mouse$NewCelltype, levels = NewCelltypesOrdered)
WT_Mouse$seurat_clusters <- factor(WT_Mouse$seurat_clusters, levels = c(0:19))
Idents(WT_Mouse) <- 'NewCelltype'
DefaultAssay(WT_Mouse) <- 'RNA'
WT_Mouse## An object of class Seurat
## 38268 features across 7101 samples within 1 assay
## Active assay: RNA (38268 features, 2500 variable features)
## 3 layers present: counts, data, scale.data
## 2 dimensional reductions calculated: pca, umapa <- DimPlot(WT_Mouse, reduction = "umap", group.by = "NewCelltype")
pbuild <- ggplot2::ggplot_build(a)
pdata <- pbuild$data[[1]]
pdata <- pdata[order(pdata$group), ]
ucols <- unique(pdata$colour)
names(ucols) <- levels(WT_Mouse$NewCelltype)
ucols## (0) E ADR Chromaffin (1) PN NOR Chromaffin (2) PN ADR Chromaffin
## "#F8766D" "#EA8331" "#D89000"
## (3) Endothelial (4) Macrophage (5) E NOR Chromaffin
## "#C09B00" "#A3A500" "#7CAE00"
## (6) E Glia|SCP (7) Cortex ZG (8) Cycling TH+
## "#39B600" "#00BB4E" "#00BF7D"
## (9) Cortex X Zone (10) E Neuroblast (11) PN Sustentacular
## "#00C1A3" "#00BFC4" "#00BAE0"
## (12) Cortex ZF (13) Macrophage (14) Cortex ZG
## "#00B0F6" "#35A2FF" "#9590FF"
## (15) Dendritic cell (16) T cell (17) NK cell
## "#C77CFF" "#E76BF3" "#FA62DB"
## (18) Mesenchymal (19) Dendritic cell
## "#FF62BC" "#FF6A98"Cell-cell-communication All Cell types
CellChat_AgedAdult <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_AgedAdult_AllCelltype_CCC.rds")
CellChat_P90 <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_3MonthAdult_AllCelltype_CCC.rds")
CellChat_P14 <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_PN_14D_AllCelltype_CCC.rds")
CellChat_P5 <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_PN_5D_AllCelltype_CCC.rds")
CellChat_E17 <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_AllCellType_SpeAge/WT_Mouse_E_17_AllCelltype_CCC.rds")CellChat_AgedAdult@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_AgedAdult@meta), "NewCelltype"])
CellChat_AgedAdult <- setIdent(CellChat_AgedAdult, ident.use = "NewCelltype")## Rename cell groups but do not change the order!## The cell group order before renaming is 1-PN_NOR Chromaffin 11-PN_Glia/Sustentacular 12-PN_Cortex ZF 15-PN_Dendritic cell 16-PN_T cell 17-PN_NK cell 18-PN_Mesenchymal 19-PN_Dendritic cell 2-PN_ADR Chromaffin 3-PN|E_Endothelial 4-PN_Macrophage 7-PN_Cortex ZG
## The cell group order after renaming is (1) PN NOR Chromaffin (11) PN Sustentacular (12) Cortex ZF (15) Dendritic cell (16) T cell (17) NK cell (18) Mesenchymal (19) Dendritic cell (2) PN ADR Chromaffin (3) Endothelial (4) Macrophage (7) Cortex ZG## Warning in setIdent(CellChat_AgedAdult, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
## These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.CellChat_P90@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_P90@meta), "NewCelltype"])
CellChat_P90 <- setIdent(CellChat_P90, ident.use = "NewCelltype")## Rename cell groups but do not change the order!## The cell group order before renaming is 1-PN_NOR Chromaffin 11-PN_Glia/Sustentacular 15-PN_Dendritic cell 16-PN_T cell 17-PN_NK cell 2-PN_ADR Chromaffin 3-PN|E_Endothelial 4-PN_Macrophage
## The cell group order after renaming is (1) PN NOR Chromaffin (11) PN Sustentacular (15) Dendritic cell (16) T cell (17) NK cell (2) PN ADR Chromaffin (3) Endothelial (4) Macrophage## Warning in setIdent(CellChat_P90, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
## These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.CellChat_P14@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_P14@meta), "NewCelltype"])
CellChat_P14 <- setIdent(CellChat_P14, ident.use = "NewCelltype")## Rename cell groups but do not change the order!## The cell group order before renaming is 1-PN_NOR Chromaffin 11-PN_Glia/Sustentacular 13-PN_Macrophage 14-PN_Cortex ZG 15-PN_Dendritic cell 16-PN_T cell 2-PN_ADR Chromaffin 3-PN|E_Endothelial 4-PN_Macrophage 7-PN_Cortex ZG 8-E_Cycling Medullary cell 9-PN_Cortex X Zone
## The cell group order after renaming is (1) PN NOR Chromaffin (11) PN Sustentacular (13) Macrophage (14) Cortex ZG (15) Dendritic cell (16) T cell (2) PN ADR Chromaffin (3) Endothelial (4) Macrophage (7) Cortex ZG (8) Cycling TH+ (9) Cortex X Zone## Warning in setIdent(CellChat_P14, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
## These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.CellChat_P5@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_P5@meta), "NewCelltype"])
CellChat_P5 <- setIdent(CellChat_P5, ident.use = "NewCelltype")## Rename cell groups but do not change the order!## The cell group order before renaming is 1-PN_NOR Chromaffin 11-PN_Glia/Sustentacular 2-PN_ADR Chromaffin 8-E_Cycling Medullary cell
## The cell group order after renaming is (1) PN NOR Chromaffin (11) PN Sustentacular (2) PN ADR Chromaffin (8) Cycling TH+## Warning in setIdent(CellChat_P5, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
## These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.CellChat_E17@meta$NewCelltype <- as.character(WT_Mouse@meta.data[rownames(CellChat_E17@meta), "NewCelltype"])
CellChat_E17 <- setIdent(CellChat_E17, ident.use = "NewCelltype")## Rename cell groups but do not change the order!## The cell group order before renaming is 0-E_ADR Chromaffin 10-E_Neuroblast 3-PN|E_Endothelial 5-E_NOR Chromaffin 6-E_Glia/SCP 7-PN_Cortex ZG 8-E_Cycling Medullary cell
## The cell group order after renaming is (0) E ADR Chromaffin (10) E Neuroblast (3) Endothelial (5) E NOR Chromaffin (6) E Glia|SCP (7) Cortex ZG (8) Cycling TH+## Warning in setIdent(CellChat_E17, ident.use = "NewCelltype"): All the calculations after `computeCommunProb` should be re-run!!
## These include but not limited to `computeCommunProbPathway`,`aggregateNet`, and `netAnalysis_computeCentrality`.CellChat_AgedAdult_WNT <- CellChat_AgedAdult
CellChat_AgedAdult_WNT@LR$LRsig <- CellChat_AgedAdult@LR$LRsig[which(CellChat_AgedAdult@LR$LRsig$ligand == "Wnt6"),]
CellChat_P90_WNT <- CellChat_P90
CellChat_P90_WNT@LR$LRsig <- CellChat_P90@LR$LRsig[which(CellChat_P90@LR$LRsig$ligand == "Wnt6"),]
CellChat_P14_WNT <- CellChat_P14
CellChat_P14_WNT@LR$LRsig <- CellChat_P14@LR$LRsig[which(CellChat_P14@LR$LRsig$ligand == "Wnt6"),]
CellChat_P5_WNT <- CellChat_P5
CellChat_P5_WNT@LR$LRsig <- CellChat_P5@LR$LRsig[which(CellChat_P5@LR$LRsig$ligand == "Wnt6"),]
CellChat_E17_WNT <- CellChat_E17
CellChat_E17_WNT@LR$LRsig <- CellChat_E17@LR$LRsig[which(CellChat_E17@LR$LRsig$ligand == "Wnt6"),]Cell-cell-communication Medullary Cell types
cellChat_Merged <- readRDS("./CellChat_CellCellCommunication/WT_Mouse_MedullaCell_SpeAge/WT_Mouse_Merged_AllAgeGroup_CCC.Rds")Figure 5 A
DimPlot(WT_Mouse, group.by='NewCelltype', reduction = 'umap', label = T, repel = T) + theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank())
Figure 5 B
Medullary Cell types
WT_Mouse_Medulla <- subset(WT_Mouse, seurat_clusters %in% c(0, 1, 2, 5, 6, 8, 10, 11))
Medulla_CellNames <- colnames(WT_Mouse_Medulla)[which((2 * WT_Mouse_Medulla@reductions$umap@cell.embeddings[,1] + WT_Mouse_Medulla@reductions$umap@cell.embeddings[,2]) > 0)]
WT_Mouse_Medulla <- subset(WT_Mouse_Medulla, cells = Medulla_CellNames)
DimPlot(WT_Mouse_Medulla, group.by='GrossAge', pt.size=0.4, reduction = 'umap', label = F, repel = T)
Figure 5 C
OutDir <- "./DEG_GO/"
WT_Mouse_SpeDEGs <- SpecificDEG(WT_Mouse, Cluster = "NewCelltype", min.logfc = 0.5, p.adj = 0.05, slot = "data", save = T, return_res = T, outdir = OutDir, Prefix = "Whole_Data_All_Cluster")Glia_DEGs <- FindMarkers(WT_Mouse, ident.1 = "(11) PN Sustentacular", ident.2 = "(6) E Glia|SCP", logfc.threshold = 0.5, group.by = "NewCelltype", only.pos = F)
saveRDS(Glia_DEGs, "Glia_Compararsion_PN_vs_E17_DEGs.rds")
head(Glia_DEGs)# Define the pattern strings
patterns <- c("^\\(11\\) PN Sustentacular", "^\\(6\\) E Glia\\|SCP")
# Function to split specified columns and filter based on the patterns
split_and_extract <- function(df, columns_to_split, patterns, columns_to_keep) {
# Initialize a copy of the dataframe to store results
df_extracted <- df[,columns_to_keep, drop = FALSE]
# Only split the specified columns
for (col in columns_to_split) {
df_extracted[, col] <- sapply(df[[col]], function(col_value) {
# Split each column element based on "| " pattern
col_split <- strsplit(col_value, "\\| ")
# Filter the values that match the patterns
filtered_values <- sapply(col_split, function(x) {
matched <- x[grepl(paste(patterns, collapse = "|"), x)]
if (length(matched) > 0) {
return(matched)
} else {
return(NA) # Return NA if no match
}
})
filtered_values <- as.numeric(gsub(paste(paste0(patterns, ": "), collapse = "|"),"",filtered_values))
return(filtered_values)
})
}
return(df_extracted)
}
# Specify the columns you want to split
columns_to_split <- c("others.pct.2", "others.pval", "others.padj", "others.log2FC") # You can change this to your column names
# Specify the columns you want to keep
columns_to_keep <- c("gene", "cluster", "pct.1") # You can change this to your column names
Genes2Show <- c("Sox10", "Wwtr1", "Tnc", "Mpz", "Lmo4", "Mef2c", "Gpr17", "Mmp2", "Ngfr", "Sox5", "Ctnnb1", "Moxd1", "Hif3a", "Atp1a2", "Olfml2a", "Col28a1", "Aspa", "Notch1", "Heyl", "Sfrp5", "Sfrp1", "Wnt6", "Gfap", "Tgfb2", "Sema3c", "Adam11", "Fam19a5", "Sparc", "Hey2", "Rgmb")
Glia_SepDEGs <- WT_Mouse_SpeDEGs[which(WT_Mouse_SpeDEGs$cluster %in% c("(11) PN Sustentacular", "(6) E Glia|SCP")),c("gene", "cluster", "pct.1", "others.pct.2", "others.pval", "others.padj", "others.log2FC")]
# Apply the function to extract the desired columns
result <- split_and_extract(Glia_SepDEGs, columns_to_split, patterns)
colnames(result) <- gsub("others.", "", colnames(result))
result$cluster2 <- ""
result$cluster2[which(result$cluster == "(6) E Glia|SCP")] <- "(11) PN Sustentacular"
result$cluster2[which(result$cluster == "(11) PN Sustentacular")] <- "(6) E Glia|SCP"
colnames(result)[2] <- "cluster1"
result <- result[,c(1,2,8,3:7)]
# View the result
# Create the volcano plot
result2 <- result
result2$log2FC[which(result2$cluster1 == "(6) E Glia|SCP")] <- -result2$log2FC[which(result2$cluster1 == "(6) E Glia|SCP")]
head(result2, 20)Glia_DEGs_M <- Glia_DEGs[which(Glia_DEGs$p_val_adj != 1),]
colnames(Glia_DEGs_M) <- c("pval", "log2FC", "pct.1", "pct.2", "padj")
Glia_DEGs_M$gene <- rownames(Glia_DEGs_M)
Glia_DEGs_M$cluster1 <- "(11) PN Sustentacular"
Glia_DEGs_M$cluster2 <- "(6) E Glia|SCP"
Glia_DEGs_M <- Glia_DEGs_M[,colnames(result2)]
Glia_DEGs_M <- Glia_DEGs_M[-which(Glia_DEGs_M$gene %in% result2$gene),]
Glia_DEGs_M$DEG_Type <- "Non-Specific"
head(Glia_DEGs_M)result2$DEG_Type <- "(6) E Glia|SCP"
result2$DEG_Type[which(result2$cluster1 == "(11) PN Sustentacular")] <- "(11) PN Sustentacular"
Glia_DEGs_Combined <- rbind(Glia_DEGs_M, result2)
Glia_DEGs_Combined$DEG_Type <- factor(Glia_DEGs_Combined$DEG_Type, levels = c("(11) PN Sustentacular", "(6) E Glia|SCP", "Non-Specific"))
# Create the volcano plot
p <- ggplot(Glia_DEGs_Combined, aes(x = log2FC, y = -log10(padj), color = DEG_Type)) + geom_point(alpha = 0.8, aes(size = ifelse(grepl("Glia-Specific", DEG_Type), 0.5, 0.2)), shape = 16) +
scale_color_manual(values = c("#01BBE1", "#3DB54A", "grey88")) + theme_minimal() + labs(title = "Volcano Plot", x = "log2 Fold Change", y = "-log10 Adjusted p-value") +
theme(legend.position = c(0.9, 0.92), plot.title = element_text(hjust = 0.5), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(color = "black"),
axis.ticks = element_line(color = "black"), legend.text = element_text(face="bold",colour = "black", size = 10), panel.border = element_rect(color = "black", fill = NA)) +
geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "black") + geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
geom_text_repel(data = result2[result2$gene %in% Genes2Show,], mapping = aes(x = log2FC, y = -log10(padj), color = DEG_Type, label = gene, size = 3),
max.overlaps = Inf, show.legend = FALSE, fontface = "bold.italic") +
theme(strip.text.x = element_text(angle=0, face="bold",size = 15), strip.text.y=element_text(angle=90,face="bold",size = 30), title = element_text(face="bold",colour = "black", size = 12)) +
theme(axis.text.x = element_text(face = "bold",colour = "black", angle = 0, hjust = 1, size = 10), axis.text.y = element_text(face="bold",colour = "black", size = 10)) +
guides(size = "none", color = guide_legend(override.aes = list(size = 4)))
p
Figure 5 D
color.heatmap <- "viridis"
color.use <- tryCatch({
RColorBrewer::brewer.pal(n = 5, name = color.heatmap)
}, error = function(e) {
(scales::viridis_pal(option = color.heatmap, direction = -1))(5)
})
MedullaCellAnnoOrder <- c("(6) E Glia|SCP", "(11) PN Sustentacular", "(10) E Neuroblast", "(8) Cycling TH+", "(0) E ADR Chromaffin", "(5) E NOR Chromaffin", "(2) PN ADR Chromaffin", "(1) PN NOR Chromaffin")
SpeDEG_Dotplot <- MyDotPlots(WT_Mouse_Medulla, result2[result2$gene %in% Genes2Show,],GeneCol = 'gene',Assay = "RNA", CelltypeCol = 'cluster1', Return = T,Draw = F, Group = "NewCelltype", CellTypeOrder = MedullaCellAnnoOrder, scale = F) +
theme(panel.background = element_blank(), plot.background = element_blank())## Loading required package: ggh4xSpeDEG_Dotplot <- SpeDEG_Dotplot + scale_color_gradientn(colors = color.use, oob = scales::squish)## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.SpeDEG_Dotplot
Figure 5 E
WT_Mouse_SCPs <- subset(WT_Mouse, NewCelltype %in% c("(6) E Glia|SCP", "(11) PN Sustentacular"))
Genes = rev(c("Hey2", "Wnt6", "Col28a1", "Notch1","Ndrg1","Sfrp5",
"Moxd1", "Sox5", "Gpr17", "Lmo4","Ngfr","Mpz",
"Serpine2","Fabp7","Erbb3","Plp1","Sox2","Sox10"))
plot_list2 <- list()
for (g in Genes) {
plot_list2[[g]] <- FeaturePlot(WT_Mouse_SCPs, features = g, order = T, alpha = 0.8, pt.size = 0.5, slot = "scale.data") +
theme(line = element_blank(), axis.text = element_blank(), axis.title = element_blank(), legend.position = "inside", legend.position.inside = c(0.1,0.7), plot.title = element_text(face = "bold.italic")) +
scale_color_gradientn(colors = c("dodgerblue", "grey90", "red"))
}## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.print(wrap_plots(plot_list2, ncol = 6) + theme(text = element_text(family = "Arial")))
Figure 5 F
OutDir <- "./DEG_GO/"
Whole_E_Glia_GO <- ComputeGO(GeneList = WT_Mouse_SpeDEGs$gene[which(WT_Mouse_SpeDEGs$cluster == "(6) E Glia|SCP")], ShowFigure = T, PlotFigure = T, SaveAsTable = T, FileName = "Whole_E_Glia_Enriched",
TopGO = 20, OrgDb = "Mouse",fromType = 'SYMBOL', toType = 'ENTREZID', Path = OutDir)## Loading required package: clusterProfiler## ## Registered S3 method overwritten by 'ggtree':
## method from
## fortify.igraph ggnetwork## clusterProfiler v4.10.1 For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
##
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141##
## Attaching package: 'clusterProfiler'## The following object is masked from 'package:igraph':
##
## simplify## The following object is masked from 'package:IRanges':
##
## slice## The following object is masked from 'package:S4Vectors':
##
## rename## The following object is masked from 'package:purrr':
##
## simplify## The following object is masked from 'package:stats':
##
## filter## Loading required package: org.Mm.eg.db## Loading required package: AnnotationDbi##
## Attaching package: 'AnnotationDbi'## The following object is masked from 'package:clusterProfiler':
##
## select## The following object is masked from 'package:dplyr':
##
## select## ## 'select()' returned 1:1 mapping between keys and columns## Warning in bitr(GeneList, fromType = fromType, toType = toType, OrgDb =
## OrgDb_Used): 6.98% of input gene IDs are fail to map...Whole_PN_Glia_GO <- ComputeGO(GeneList = WT_Mouse_SpeDEGs$gene[which(WT_Mouse_SpeDEGs$cluster == "(11) PN Sustentacular")], ShowFigure = T, PlotFigure = T, SaveAsTable = T, FileName = "Whole_PN_Glia_Enriched",
TopGO = 20, OrgDb = "Mouse",fromType = 'SYMBOL', toType = 'ENTREZID', Path = OutDir)## 'select()' returned 1:1 mapping between keys and columns## Warning in bitr(GeneList, fromType = fromType, toType = toType, OrgDb =
## OrgDb_Used): 2.22% of input gene IDs are fail to map...E17_Glia_GO <- read.table(paste0(OutDir,"Whole_E_Glia_Enriched_GO_Sig.tsv"), header = T, sep = "\t")
PN_Glia_GO <- read.table(paste0(OutDir,"Whole_PN_Glia_Enriched_GO_Sig.tsv"), header = T, sep = "\t")
PN_Glia_GO_Spe <- PN_Glia_GO[-which(PN_Glia_GO$ID %in% E17_Glia_GO$ID),]
P_PNGlia_PathwayGO <- Bubble_Plot(pathway = PN_Glia_GO_Spe[grep("signaling", PN_Glia_GO_Spe$Description),][1:25,], "PN Glia Cell All Pathways")
P_PNGlia_PathwayGO + theme(text = element_text(family = "Arial"))
Figure 5 G
ucols_subset <- ucols[levels(CellChat_P90@idents)]
netVisual_individual(CellChat_P90, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)
## [[1]]Figure 5 H
ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
netVisual_individual(CellChat_AgedAdult, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)
## [[1]]Figure 5 I
gg <- netVisual_bubble(object = cellChat_Merged, sources.use = c(1:3), targets.use = c(1:3), comparison = c(1, 2, 3, 4, 5), angle.x = 90, signaling = c("WNT", "NOTCH"), return.data = T)## Comparing communications on a merged objectSelectedCommunicationDF <- gg$communication
SelectedCommunicationDF_Modified <- SelectedCommunicationDF[grep("Wnt",SelectedCommunicationDF$interaction_name_2, invert = T),]
WNT_df <- SelectedCommunicationDF[grep("Wnt", SelectedCommunicationDF$interaction_name_2),]
WNT_df_Modified <- WNT_df %>% na.omit() %>% group_by(source.target, ligand) %>% reframe(
prob = sum(prob, na.rm = TRUE),
prob.original = sum(prob.original, na.rm = TRUE),
pval = round(mean(pval, na.rm = TRUE)),
interaction_name = paste(ligand, "- All Activated Receptors"),
interaction_name_2 = paste(ligand, "- All Activated Receptors"),
receptor = "All Activated Receptors",
across(-c(prob, prob.original, pval, interaction_name, interaction_name_2), dplyr::first)) %>% ungroup()
WNT_df_Modified <- WNT_df_Modified[,colnames(SelectedCommunicationDF)]
SelectedCommunicationDF_Modified <- rbind(SelectedCommunicationDF_Modified, WNT_df_Modified)g <- Figure5I_Quick(CellChatObj = cellChat_Merged, CommunicationDF = SelectedCommunicationDF_Modified)
g + theme(text = element_text(family = "Arial"))
Figure 5 J
Notch_Genes <- read.table("./Notch_Genes.csv", header = T, sep=",")
Notch_Genes$Gene <- DescTools::StrCap(Notch_Genes$Gene, method = "title")## Registered S3 methods overwritten by 'proxy':
## method from
## print.registry_field registry
## print.registry_entry registryMyDotPlots(WT_Mouse_Medulla, Notch_Genes[which(Notch_Genes$Fine_Type != "Diff"),],GeneCol = 'Gene', Assay = "RNA", CelltypeCol = 'Fine_Type',Return = T,Draw = F, Group = "NewCelltype", CellTypeOrder = MedullaCellAnnoOrder, scale = F) +
theme(panel.background = element_blank(), plot.background = element_blank()) + scale_color_gradientn(colors = color.use, oob = scales::squish)## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
Sup Figure 5 A - C
tmp1 <- DimPlot(WT_Mouse, group.by=c("samples"), reduction = 'umap', label = F) + theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank()) + guides(color = guide_legend(ncol = 1))
tmp2 <- DimPlot(WT_Mouse, group.by=c("DetailedAge"), reduction = 'umap', label = F) + theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank())
tmp3 <- DimPlot(WT_Mouse, group.by=c("seurat_clusters"), reduction = 'umap', label = T, repel = T) + theme(axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank()) + NoLegend()
tmp1 + tmp2 + tmp3
Sup Figure 5 D
MyMarkers <- read.table("./Mouse_AdrenalGland_Markers.txt", header = T, sep="\t")
MyMarkers$Markers <- DescTools::StrCap(MyMarkers$Markers, method = "title")
MyDotPlots(WT_Mouse, MyMarkers,GeneCol = 'Markers',Assay = "RNA", CelltypeCol = 'Celltype',Return = T,Draw = F, ClusterIdents = F, scale = F,
CellTypeOrder = c("6","11","10","8","0","2","1","5","18","7","14","12","9","3","4","13","16","15","19","17")) +
theme(panel.background = element_blank(), plot.background = element_blank(), axis.text.x = element_text(vjust = 0.5, face = "italic"), axis.text.y = element_text(face = "plain"))
Sup Figure 5 E
ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
netAnalysis_signalingRole_network(CellChat_AgedAdult, signaling = "NOTCH", width = 10, height = 3, font.size = 10, color.use = ucols_subset)
Sup Figure 5 F
ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
netAnalysis_signalingRole_network(CellChat_AgedAdult, signaling = "WNT", width = 10, height = 3, font.size = 10, color.use = ucols_subset)
Sup Figure 5 G
ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
netAnalysis_signalingRole_network(CellChat_AgedAdult, signaling = "BMP", width = 10, height = 3, font.size = 10, color.use = ucols_subset)
Sup Figure 5 H
layout(matrix(1:5, nrow = 1, byrow = TRUE))
ucols_subset <- ucols[levels(CellChat_E17@idents)]
p1 <- netVisual_individual(CellChat_E17, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)
ucols_subset <- ucols[levels(CellChat_P5@idents)]
p2 <- netVisual_individual(CellChat_P5, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)
ucols_subset <- ucols[levels(CellChat_P14@idents)]
p3 <- netVisual_individual(CellChat_P14, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)
ucols_subset <- ucols[levels(CellChat_P90@idents)]
p4 <- netVisual_individual(CellChat_P90, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)
ucols_subset <- ucols[levels(CellChat_AgedAdult@idents)]
p5 <- netVisual_individual(CellChat_AgedAdult, signaling = "NOTCH", pairLR.use = "DLK1_NOTCH1", layout = "circle", color.use = ucols_subset)
Sup Figure 5 I
layout(matrix(1:4, nrow = 1, byrow = TRUE))
ucols_subset <- ucols[levels(CellChat_P5_WNT@idents)]
netVisual_aggregate(CellChat_P5_WNT, signaling = "WNT", layout = "circle", color.use = ucols_subset)
ucols_subset <- ucols[levels(CellChat_P14_WNT@idents)]
netVisual_aggregate(CellChat_P14_WNT, signaling = "WNT", layout = "circle", color.use = ucols_subset)
ucols_subset <- ucols[levels(CellChat_P90_WNT@idents)]
netVisual_aggregate(CellChat_P90_WNT, signaling = "WNT", layout = "circle", color.use = ucols_subset)
ucols_subset <- ucols[levels(CellChat_AgedAdult_WNT@idents)]
netVisual_aggregate(CellChat_AgedAdult_WNT, signaling = "WNT", layout = "circle", color.use = ucols_subset) 
Session Information
sessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: SUSE Linux Enterprise Server 15 SP5
##
## Matrix products: default
## BLAS/LAPACK: /cfs/klemming/projects/supr/sllstore2017016/joey.d/Tools/anaconda3/envs/R432/lib/libopenblasp-r0.3.27.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Stockholm
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] org.Mm.eg.db_3.18.0 AnnotationDbi_1.64.1
## [3] clusterProfiler_4.10.1 ggh4x_0.2.8
## [5] showtext_0.9-7 showtextdb_3.0
## [7] sysfonts_0.8.9 CellChat_1.6.1
## [9] ggpubr_0.6.0 reshape2_1.4.4
## [11] mclust_6.1.1 cluster_2.1.6
## [13] cowplot_1.1.3 ggthemes_5.1.0
## [15] ComplexHeatmap_2.18.0 gridGraphics_0.5-1
## [17] gridExtra_2.3 extrafont_0.19
## [19] ggrepel_0.9.5 patchwork_1.2.0
## [21] igraph_2.0.3 network_1.19.0
## [23] viridis_0.6.5 viridisLite_0.4.2
## [25] SCopeLoomR_0.13.0 SingleCellExperiment_1.24.0
## [27] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [29] GenomicRanges_1.54.1 GenomeInfoDb_1.38.8
## [31] IRanges_2.36.0 S4Vectors_0.40.2
## [33] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [35] matrixStats_1.3.0 magrittr_2.0.3
## [37] lubridate_1.9.3 forcats_1.0.0
## [39] stringr_1.5.1 purrr_1.0.2
## [41] readr_2.1.5 tidyr_1.3.1
## [43] tibble_3.2.1 ggplot2_3.5.1
## [45] tidyverse_2.0.0 Seurat_5.0.0
## [47] SeuratObject_5.0.2 sp_2.1-1
## [49] dplyr_1.1.4 openxlsx_4.2.7
##
## loaded via a namespace (and not attached):
## [1] gld_2.6.6 goftest_1.2-3 Biostrings_2.70.3
## [4] vctrs_0.6.5 spatstat.random_3.2-3 proxy_0.4-27
## [7] digest_0.6.35 png_0.1-8 shape_1.4.6.1
## [10] Exact_3.3 registry_0.5-1 deldir_2.0-4
## [13] parallelly_1.37.1 magick_2.8.3 MASS_7.3-60
## [16] httpuv_1.6.15 foreach_1.5.2 qvalue_2.34.0
## [19] withr_3.0.0 xfun_0.49 ggfun_0.1.5
## [22] survival_3.7-0 memoise_2.0.1 gson_0.1.0
## [25] systemfonts_1.1.0 ragg_1.3.2 tidytree_0.4.6
## [28] zoo_1.8-12 GlobalOptions_0.1.2 pbapply_1.7-2
## [31] KEGGREST_1.42.0 promises_1.3.0 httr_1.4.7
## [34] rstatix_0.7.2 globals_0.16.3 fitdistrplus_1.1-11
## [37] rstudioapi_0.16.0 miniUI_0.1.1.1 generics_0.1.3
## [40] DOSE_3.28.2 ggalluvial_0.12.5 zlibbioc_1.48.2
## [43] ggraph_2.2.1 polyclip_1.10-6 GenomeInfoDbData_1.2.11
## [46] SparseArray_1.2.4 xtable_1.8-4 doParallel_1.0.17
## [49] evaluate_0.24.0 S4Arrays_1.2.1 hms_1.1.3
## [52] bookdown_0.41 irlba_2.3.5.1 colorspace_2.1-0
## [55] hdf5r_1.3.11 ggnetwork_0.5.13 ROCR_1.0-11
## [58] readxl_1.4.3 reticulate_1.38.0 spatstat.data_3.0-4
## [61] lmtest_0.9-40 ggtree_3.10.1 later_1.3.2
## [64] lattice_0.22-6 spatstat.geom_3.2-9 NMF_0.28
## [67] future.apply_1.11.2 shadowtext_0.1.4 scattermore_1.2
## [70] RcppAnnoy_0.0.22 class_7.3-22 pillar_1.9.0
## [73] nlme_3.1-165 iterators_1.0.14 sna_2.8
## [76] gridBase_0.4-7 compiler_4.3.2 RSpectra_0.16-1
## [79] stringi_1.8.4 DescTools_0.99.58 tensor_1.5
## [82] plyr_1.8.9 crayon_1.5.2 abind_1.4-5
## [85] haven_2.5.4 graphlayouts_1.2.0 bit_4.0.5
## [88] rootSolve_1.8.2.4 fastmatch_1.1-4 textshaping_0.3.7
## [91] codetools_0.2-20 bslib_0.7.0 e1071_1.7-14
## [94] lmom_3.2 GetoptLong_1.0.5 plotly_4.10.4
## [97] mime_0.12 splines_4.3.2 circlize_0.4.16
## [100] Rcpp_1.0.12 fastDummies_1.7.3 HDO.db_0.99.1
## [103] cellranger_1.1.0 Rttf2pt1_1.3.12 knitr_1.47
## [106] blob_1.2.4 utf8_1.2.4 clue_0.3-65
## [109] fs_1.6.4 listenv_0.9.1 expm_1.0-0
## [112] ggsignif_0.6.4 ggplotify_0.1.2 Matrix_1.6-5
## [115] statmod_1.5.0 tzdb_0.4.0 svglite_2.1.3
## [118] tweenr_2.0.3 pkgconfig_2.0.3 tools_4.3.2
## [121] cachem_1.1.0 RSQLite_2.3.7 DBI_1.2.3
## [124] fastmap_1.2.0 rmarkdown_2.27 scales_1.3.0
## [127] ica_1.0-3 broom_1.0.6 sass_0.4.9
## [130] coda_0.19-4.1 FNN_1.1.4 BiocManager_1.30.23
## [133] dotCall64_1.1-1 carData_3.0-5 RANN_2.6.1
## [136] farver_2.1.2 scatterpie_0.2.4 tidygraph_1.3.1
## [139] yaml_2.3.8 cli_3.6.2 leiden_0.4.3.1
## [142] lifecycle_1.0.4 uwot_0.2.2 mvtnorm_1.2-5
## [145] presto_1.0.0 backports_1.5.0 BiocParallel_1.36.0
## [148] timechange_0.3.0 gtable_0.3.5 rjson_0.2.21
## [151] ggridges_0.5.6 progressr_0.14.0 ape_5.8
## [154] parallel_4.3.2 limma_3.58.1 jsonlite_1.8.8
## [157] RcppHNSW_0.6.0 bitops_1.0-7 bit64_4.0.5
## [160] Rtsne_0.17 yulab.utils_0.1.8 spatstat.utils_3.0-5
## [163] BiocNeighbors_1.20.2 zip_2.3.1 jquerylib_0.1.4
## [166] highr_0.11 GOSemSim_2.28.1 lazyeval_0.2.2
## [169] shiny_1.8.1.1 htmltools_0.5.8.1 enrichplot_1.22.0
## [172] GO.db_3.18.0 sctransform_0.4.1 glue_1.7.0
## [175] spam_2.10-0 XVector_0.42.0 RCurl_1.98-1.14
## [178] treeio_1.26.0 boot_1.3-30 extrafontdb_1.0
## [181] R6_2.5.1 labeling_0.4.3 rngtools_1.5.2
## [184] aplot_0.2.3 statnet.common_4.10.0 DelayedArray_0.28.0
## [187] tidyselect_1.2.1 ggforce_0.4.2 car_3.1-2
## [190] future_1.33.2 munsell_0.5.1 KernSmooth_2.23-24
## [193] data.table_1.15.4 htmlwidgets_1.6.4 fgsea_1.28.0
## [196] RColorBrewer_1.1-3 rlang_1.1.4 spatstat.sparse_3.0-3
## [199] spatstat.explore_3.2-7 Cairo_1.6-2 fansi_1.0.6